

%% load estimation results
load('parameter_estimates_direction_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z specif


%% Load everything of interest
load('sim_eq_direction_round2.mat')
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% standard errors via Hessian
InvSigmaH = HessPhat(thetaStar,1,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
seH = sqrt(diag(inv(InvSigmaH)));
resultsH = printResults([thetaStar,seH,2*normcdf(abs(thetaStar./seH),'upper')],reb,mpow,specif);
disp(' ')
disp('for Tables F1, F3-F4:')
disp(resultsH)


%%
clearvars -except thetaStar InvSigmaH seH resultsH
save('inference_direction.mat')



